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Abstract 

We describe an elementary algorithm to build convex inner approximations of 
nonconvex sets. Both input and output sets are basic semialgebraic sets given as 
lists of defining multivariate polynomials. Even though no optimality guarantees 
can be given (e.g. in terms of volume maximization for bounded sets), the algorithm 
is designed to preserve convex boundaries as much as possible, while removing re- 
gions with concave boundaries. In particular, the algorithm leaves invariant a given 
convex set. The algorithm is based on Gloptipoly 3, a public-domain Matlab pack- 
age solving nonconvex polynomial optimization problems with the help of convex 
semidefinite programming (optimization over linear matrix inequalities, or LMIs). 
We illustrate how the algorithm can be used to design fixed-order controllers for 
linear systems, following a polynomial approach. 
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1 Introduction 



The set of controllers stabilizing a linear system is generally nonconvex in the parameter 
space, and this is an essential difficulty faced by numerical algorithms of computer-aided 
control system design, see e.g. [4] and references therein. It follows from the derivation 
of the Routh-Hurwitz stability criterion (or its discrete-time counterpart) that the set 
of stabilizing controllers is real basic semialgebraic, i.e. it is the intersection of sublevel 
sets of given multivariate polynomials. A convex inner approximation of this nonconvex 
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semialgebraic stability region was obtained in [I] in the form of linear matrix inequali- 
ties (LMI) obtained from univariate polynomial positivity conditions, see also [9] . Convex 
polytopic inner approximations were also obtained in [TB] , for discrete-time stability, using 
reflection coefficients. Convex inner approximations make it possible to design stabiliz- 
ing controllers with the help of convex optimization techniques, at the price of loosing 
optimality w.r.t. closed-loop performance criteria (if 2 norm, norm or alike). 

Generally speaking, the technical literature abounds of convex outer approximations of 
nonconvex semialgebraic sets. In particular, such approximations form the basis of many 
branch-and-bound global optimization algorithms [15) . By construction, Lasserre's hier- 
archy of LMI relaxations for polynomial programming is a sequence of embedded convex 
outer approximations which are semidefinite representable, i.e. which are obtained by 
projecting affme sections of the convex cone of positive semidefinite matrices, at the price 
of introducing lifting variables [6] . 

After some literature search, we could not locate any systematic constructive procedure 
to generate convex inner approximations of nonconvex semialgebraic sets, contrasting 
sharply with the many convex outer approximations mentioned above. In the context of 
fixed-order controller design, inner approximations correspond to a guarantee of stability, 
at the price of loosing optimality. No such stability guarantee can be ensured with outer 
approximations. 

The main contribution of this paper is therefore an elementary algorithm, readily imple- 
mentable in Matlab, that generates convex inner approximations of nonconvex sets. Both 
input and output sets are basic semialgebraic sets given as lists of defining multivariate 
polynomials. Even though no optimality guarantees can be given in terms of volume 
maximization for bounded sets, the algorithm is designed to preserve convex boundaries 
as much as possible, while removing regions with concave boundaries. In particular, the 
algorithm leaves invariant a given convex set. The algorithm is based on Gloptipoly 3, 
a public-domain Matlab package solving nonconvex polynomial optimization problems 
with the help of convex LMIs [7J. Even though the algorithm can be useful on its own, 
e.g. for testing convexity of semialgebraic sets, we illustrate how it can be used to design 
fixed-order controllers for linear systems, following a polynomial approach. 

2 Convex inner approximation 

Given a basic closed semialgebraic set 



where pi are multivariate polynomials, we are interested in computing another basic closed 
semialgebraic set 




(1) 




(2) 



which is a valid convex inner approximation of S, in the sense that 



S c S. 



7 



Ideally, we would like to find the tightest possible approximation, in the sense that the 
complement set S\S = {x G S : x S} is as small as possible. Mathematically we may 
formulate the problem as the volume minimization problem 




but since set S is not necessarily bounded we should make sure that this integral makes 
sense. Moreover, computing the volume of a given semialgebraic set is a difficult task 
in general [8], so we expect that optimizing such a quantity is as much as difficult. In 
practice, in this paper, we will content ourselves of an inner approximation that removes 
the nonconvex parts of the boundary and keeps the convex parts as much as possible. 



3 Detecting nonconvexity 



Before describing the method, let us recall some basics definitions on polynomials and 
differential geometry. Let x G M n h- > Pi(x) G M[x] be a multivariate polynomial of total 
degree d. Let 

dpi(x) 



9i{x) 



dxj 



G R n tol 



be its gradient vector and 



HAx) 



J Jj=l...n 



d 2 pi(x) 

dxjdx kJjl , n 



G R" xrt [d 



its (symmetric) Hessian polynomial matrix. Define the optimization problem 

qi = mm Xiy y T Hi(x,y)y 
s.t. Pi(x) = 

Pj (x) < 0, j = 1 . . .m, j ^ i 
y T 9i(x) = 

y T y = i 

with global minimizers {x l . . . x kl } and {y 1 . . . y kt }. 

Let us make the following nondegeneracy assumption on defining polynomials Pi(x): 



(3) 



Assumption 1 There is no point x such that Pi(x) and gi(x) vanish simultaneously while 
satisfying Pj(x) < for j — 1, . . . , m, j ^ i. 

Since the polynomial system Pi(x) = 0, gi(x) = 0, involves n + 1 equations for n un- 
knowns, Assumption [1] is satisfied generically. In other words, in the Euclidean space of 
coefficients of polynomials Pi(x), instances violating Assumption [1] belong to a variety of 
Lebesgue measure zero, and an arbitrarily small perturbation on the coefficients generates 
a perturbed set S e satisfying Assumption [TJ 
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Theorem 1 Under Assumption^ polynomial level set (CEJ) is convex if and only if qi > 
for all i = 1, . . . , m. 

Proof: The boundary of set S consists of points x such that Pi(x) = for some i, and 
Pji x ) — f° r 3 7^ i- m the neighborhood of such a point, consider the Taylor series 

Pi(x + y)= Pi{x) + y T gi(x) + y T Hi{x)y + 0(y 3 ) (4) 

where 0(y 3 ) denotes terms of degree 3 or higher in entries of vector y, the local coordinates. 
By Assumption [TJ the gradient gi(x) does not vanish along the boundary, and hence 
convexity of the boundary is inferred from the quadratic term in expression More 
specifically, when y T gi(x) = 0, vector y belongs to the hyperplane tangent to S at point 
x. Let V be a matrix spanning this linear subspace of dimension n — 1 so that y = 
Vy for some y. The quadratic form y T Hi(x)y = y T V T HiVy can be diagonalised with 
the congruence transformation y = Uy (Schur decomposition), and hence y T Hi(x)y = 
y T U T V T HiVUy T = YHiZi hi{x)yf. The eigenvalues hi(x), i — 1, . . . , n — 1 are reciprocals 
of the principal curvatures of the surface. Problem ([3]) then amounts to finding the 
minimum curvature, which is non-negative when the surface is locally convex around x.U 

In the case of three-dimensional surfaces (n = 3), the ideas of tangent plane, local coor- 
dinates and principal curvatures used in the proof of Theorem [T] are standard notions of 
differential geometry, see e.g. Section 3.3. in [2] for connections between principal cur- 
vatures and eigenvalues of the local Hessian form (called the second fundamental form, 
once suitably normalized). 




Figure 1: Hyperboloid of one sheet (white), with tangent plane (gray) at the origin, a 
saddle point with a tangent convex parabola (thick black) and a tangent concave hyperbola 
(thick black). 

As an example illustrating the proof of Theorem [H consider the hyperboloid of one sheet 
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S = {x e M 3 : 



Pi(x) = x\ 



x 2 ~~ x 3 < 0} with gradient and Hessian 



2xi 
-2x 2 
-1 




2 













2 









At the origin x = 0, the tangent plane is T = {y e M 3 : y 3 = 0} and = 2?/^ — 2y\ 

is a bivariate quadratic form with eigenvalues 2 and —2, corresponding respectively to 
the convex parabola {x : x| + x 3 = 0} (positive curvature) and concave hyperbola 
{x : x\ — X3 = 0} (negative curvature), see Figure [TJ 

Theorem [1] can be exploited in an algorithmic way to generate a convex inner approxima- 
tion of a semialgebraic set. 

Algorithm 1 (Convex inner approximation) 

Input: Polynomials pi, i — 1 . . . m defining set S as in (Op. Small nonnegative scalar e. 
Output: Polynomials pi, i = 1 . . . rh defining set S as in (TJ^. 
Step 1 : Let i — 1 . 

Step 2: If deg Pi < 1 then go to Step 5. 

Step 3: If pi(x) G S, solve optimization problem (Tj|) for optimum qi and minimizers 
{x 1 . . . x k }. If pi(x) ^ S, go to Step 5. 

Step 4: If q { < 0, then select one of the minimizers x^ , j = 1 . . . k i} let p m+ \ = gi(x^)(x — 
x J ) + e. Then let m = m + 1, and go to step 3. 

Step 5: Let % — i + 1 . If i < m then go to Step 2. 

Step 6: Return pi = p i7 i = 1, . . . m. 

The idea behind the algorithm is as follows. At Step 3, by solving the polynomial opti- 
mization problem of Theorem [1] we identify a point of minimal curvature along algebraic 
varieties defining the boundary of S. If the minimal curvature is negative, then we sep- 
arate the point from the set with a gradient hyperplane, and we iterate on the resulting 
semialgebraic set. At the end, we obtain a valid inner approximation. 

Note that Step 2 checks if the boundary is affine, in which case the minimum curvature 
is zero and there is no optimization problem to be solved. 

The key parameter of the algorithm is the small positive scalar e used at Step 4 for 
separating strictly a point of minimal curvature, so that the algorithm does not identify 
it again at the next iteration. Moreover, in Step 4, one must elect arbitrarily a minimizer. 
We will discuss this issue later in this paper. 

Finally, as pointed out to us by a referee, the ordering of the sequence of input polynomials 
Pi has an impact on the sequence of output polynomials pj, and especially on the size of 
the convex inner approximation S. However, it seems very difficult to design a priori an 
optimal ordering policy. 



4 Matlab code and geometric examples 



At each step of Algorithm [T] we have to solve a potentially nonconvex polynomial optimiza- 
tion problem. For that purpose, we use Gloptipoly 3, a public-domain Matlab package [7]. 
The methodology consists in building and solving a hierarchy of embedded linear matrix 
inequality (LMI) relaxations of the polynomial optimization problem, see the survey [12] . 
The LMI problems are solved numerically with the help of any semidefinite programming 
solver (by default Gloptipoly 3 uses SeDuMi). Under the assumption that our original 
semi- algebraic set is compact, the sequence of minimizers obtained by solving the LMI 
relaxations is ensured to converge mononotically to the global minimum. Under the addi- 
tional assumption that the global optima live on a zero-dimensional variety (i.e. there is 
a finite number of them), Gloptipoly 3 eventually extracts some of them (not necessarily 
all, but at least one) using numerical linear algebra. The LMI problems in the hierarchy 
have a growing number of variables and constraints, and the main issue is that we cannot 
predict in advance how large has to be the LMI problem to guarantee global optimality. 
In practice however we observe that it is not necessary to go very deep in the hierarchy 
to have a numerical certificate of global optimality. 

4.1 Hyperbola 

Let us first with the elementary example of an unbounded nonconvex hyperbolic region 
S = {x G IR 2 : pi(x) < 0} with p\{x) = — 1 + X\X2, for which optimization problem ((2]) 
reads 

min 2y x y 2 

s.t. x 2 yi + x x V2 = 

— 1 + X\X 2 = 

yf + y! = i. 

Necessary optimality conditions yield immediately k\ — 2 global minimizers x 1 = ^(1,1), 

y 1 = 2 (-*■' — 1) an d x<2 = "^r( — 1> — 1); y 2 = "^r( — 1> 1)> an d hence two additional (normal- 
ized) affine constraints p 2 (x) = —2 + x\ + X2 and p%{x) = —2 — x\ — X2 defining the slab 
S = {x : Pi(x) < 0, % = 1, 2, 3} = {x : —2 < x% + X2 < 2} which is indeed a valid inner 
approximation of S. 

4.2 Egg quart ic 

Now we show that Algorithm [1] can be used to detect convexity of a semialgebraic set. 
Consider the smooth quartic sublevel set S = {x e M 2 : pi(x) = x\ + x\ + x\ + X2 < 0} 
represented on Figure [2j Assumption [1] is ensured since the gradient gi(x) = [2xx(xf + 
2) 4^2 + 1] cannot vanish for real x. 

A Matlab implementation of the first steps of the algorithm can be easily written using 
Gloptipoly 3: 

% problem data 



- 



-0.2 - 



-0.4 - 



-0.6 - 



-0.8 - 



-1 - 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

x 1 

Figure 2: Convex smooth quartic. 

mpol x y 2 

pi = x(l)"4+x(2)~4+x(l)"2+x(2) ; 
gl = dif f (pi ,x) ; 7, gradient 
HI = diff(gl.x); 7. Hessian 
7« LMI relaxation order 
order = 3; 

°/ build LMI relaxation 

P = msdp(min(y'*Hl*y) , pl==0, ... 

gl*y==0, y J *y==l , order); 
°/ solve LMI relaxation 
[status, obj] = msol(P) 

Notice that we solve the LMI relaxation of order 3 (e.g. moments of degree 6) of problem 
In GloptiPoly, an LMI relaxation is solved with the command msol which returns 
two output arguments: status and obj. Argument status can take the following values: 

• -1 if the LMI relaxation is infeasible or could not be solved for numerical reasons; 

• if the LMI relaxation could be solved but it is impossible to detect global optimality 
and to extract global optimizers, in which case obj is a lower (resp. upper) bound 
on the global minimum (resp. maximum) of the original optimization problem; 

• +1 if the LMI relaxation could be solved, global optimality is certified and global 
minimizers are extracted, in which case obj is the global optimum of the original 
optimization problem. 
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Running the above script, Gloptipoly returns obj = 2.0000 and status = 1, certifying 
that the minimal curvature is strictly positive, and hence that the polynomial sublevel 
set is convex. 

Note that in this simple case, convexity of set S follows directly from positive semidefi- 
niteness of the Hessian H\(x) = diag (12x 2 + 2, 12x?,), yet Algorithm [1] can systematically 
detect convexity in more complicated cases. 

4.3 Waterdrop quartic 

Consider the quartic S = {x G 1R 2 : Pi(x) = x\ + x\ + x\ + x\ < 0} which has a singular 
point at the origin, hence violating Assumption [TJ 

Applying Algorithm [H the LMI relaxation of order 4 (moments of degree 8) yields a glob- 
ally minimal curvature of —0.094159 achieved at the 2 points x 1 = (—0.048892, —0.14076) 
and x 2 = (0.048896, —0.14076). With the two additional affme constraints Pk{x) = 
gi(x k )(x — x k ) < 0, k = 2,3, the resulting set S has a globally minimal curvature of 1 
certified at the LMI relaxation of order 4, and therefore it is a valid convex inner approx- 
imation of S, see Figure [31 
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x 1 

Figure 3: Nonconvex waterdrop quartic (light gray) and its convex inner approximation 
(dark gray) obtained by adding affme constraints at two points x 1 and x 2 of minimal 
curvature. 

This example illustrates that Algorithm [T] can work even when Assumption [T] is violated. 
Here the singularity is removed by the additional affine constraints. This example also 
shows that symmetry of the problem can be exploited, since two global minimizers are 
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found (distinct points with the same minimal curvature) to remove two nonconvex parts 
of the boundary simultaneously. 

4.4 Singular quart ic 

Consider the quartic S = {x 6 1R 2 : P\(x) = x\ + x\ + x\ < 0} which has a singular point 
at the origin, hence violating Assumption [TJ 

Running Algorithm (TJ we obtain the following sequence of bounds on the minimum cur- 
vature, for increasing LMI relaxation orders: 



order 


2 


3 


4 


5 


obj 


-7.5000 • 10- 1 


-7.7502 ■ 10- 2 


-8.5855 ■ 10~ 3 


-4.9525 ■ lO" 3 



GloptiPoly is not able to certify global optimality, so we can only speculate that the global 
minimum is zero and hence that set S is convex, see Figure HI We may say that set S is 
numerically convex. 



o - 



-0.2 - 



-0.4 - 



-0.6 - 



-0.8 - 



-1 - 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

x 1 

Figure 4: Numerically convex singular quartic. 

Indeed if we strenghten the constraint pi(x) < into pi(x)+e < for a small positive e, say 
10~ 3 , then GloptiPoly 3 certifies global optimality and convexity with obj = -4 . 0627e-7 
at the 4th LMI relaxation. On the other hand, if we relax the constraint into pi(x) +e < 
with a negative e = — 10~ 3 , then GloptiPoly 3 certifies global optimality and nonconvexity 
with obj = -0.22313 at the 4th LMI relaxation. We can conclude that the optimum of 
problem E] is sensitive, or ill-conditioned, with respect to the problem data, the coefficients 
of pi(x). The reason behind this ill-conditioning is the singularity of S at the origin, see 
Figure |5] which represents the effect of perturbing the constraint p\ (x) < around the 
singularity. 



0.15 i 1 1 1 1 1 1 0.15 




-0.2 - -0.2 - 



-0.25 1 1 1 ' ' 1 1 -0.25 1 ' ' ' 1 1 

-0.2 -0.1 0.1 0.2 -0.2 -0.1 0.1 0.2 

x 1 x 1 

Figure 5: Perturbed quartic pi(x) + e < (bold line) can be convex (e = 1CT 3 ) or 
nonconvex (e = — 10~ 3 ) near singularity of original quartic level set pi(x) = (light line). 

5 Control applications 

In this section we focus on control applications of Algorithm [I] which is used to generate 
convex inner approximation of stability regions in the parameter space. 

5.1 Third-order discrete-time stability region 

Algorithm [1] can lend insight into the (nonconvex) geometry of the stability region. Con- 
sider the simplest non-trivial case of a third-order discrete-time polynomial X\ + x 2 z + 
X3Z 2 + z 3 which is stable (roots within the open unit disk) if and only if parameter 
x = (xx,X2,xs) lies within the interior of compact region S = {x 6 M 3 : p\(x) = 
—xi - x 2 - x 3 - 1 < 0, p 2 (x) = x\ - x 2 + x 3 - 1 < 0, pz{x) —x\ — xxx?, + x 2 - 1 < 0}. 
Stability region S is nonconvex, delimited by two planes Pi(x) = 0, p 2 (x) = and a 
hyperbolic paraboloid ^3(0;) = see e.g. [U Example 11.4]. 

Optimization problem ()3]) corresponding to convexity check of the hyperbolic paraboloid 



in 



reads as follows: 

min -2y\ + 2j/ij/ 3 

S.t. x\ — X1X3 + X 2 — 1 = 

-£l - X 2 - X 3 - 1 < , s 

^1 — ^2 + ^3 — 1 < 

(2a;i - £3)2/1 + y 2 + 2*2/3 = 

yf + yi + vi = i- 

The objective function and the last constraint depend only on and necessary optimality 
conditions obtained by differentiating the Lagrangian —2y\ + 2yiy 3 + t(yf + y\ + y\ — 1) 
with respect to y yield the symmetric pencil equation 



-4 + 2t 





2 











2t 







1/2 


= 


2 





2t 




1/3 





From the determinant of the above 3-by-3 matrix, equal to t(t 2 — 2t — 1), we conclude 
that multiplier t can be equal to 1 — \/2,0 or 1 + a/2. The choice t — implies ?/i = 
0, 2/2 = 1; 2/3 = which is inconsistent with the last but one constraint in (J5]). The choice 
t = 1 - a/2 yields j/ x = ±(1 + a/2)«, y 2 = 0, y 3 = ±a with a = l/y/4- 2y/2 and the 
objective function — 2y\ + 2yiy 3 = — 1 + a/2. The choice t = 1 + a/2 yields ?/i = ±a, 
1/2 = 0, ?/3 = ±(— 1 — v^)^ and the objective function —1 — y/2, a negative minimum 
curvature. Therefore region S is indeed nonconvex. 

From the remaining constraints in fl5]), we conclude that the minimal curvature points x 
can be found along the portion of parabola y/2xf — X2 + 1 = included in the half-planes 
(2 + ^/2)£i+X2 + l > and — (2 + -\/2)2;i + S2 + l > 0. Any plane tangent to the hyperbolic 
paraboloid pz(x) = at a point along the parabola \/2x\ — x 2 + 1 = can be used to 
generate a valid inner approximation of the stability region. For example, with the choice 
x l = (0, 1, 0), we generate the gradient half-plane Pi(x) = g 3 (x r )(x — x 1 ) = x 2 — 1 < 0. 

More generally, for discrete-time polynomials of degree n > 3, stability region S is the 
image of the box B = [—1,1]™ (of so-called reflection coefficients) though a multiaffine 
mapping, see e.g. pSj and references therein. The boundary of S consists of ruled 
surfaces, and the convex hull of S is generated by the images of the vertices of B through 
the multiaffine mapping. It would be interesting to investigate whether this particular 
geometry can be exploited to generate systematically a convex inner approximation of 
maximum volume of the stability region S. 

5.2 Fixed-order controller design 

Consider the open-loop discrete-time system (— 2z 2 + l)/(z 3 + z 2 + a), parametrized by 
a G R, in negative closed-loop configuration with the controller [x,\z + x<i)j{z + 1). The 
characteristic polynomial is equal to q(z) = Y^t=o Qk% k = z 4 + 2(1 — xi)z 3 + (1 — 2x 2 )z 2 + 
(a + X\)z + a + x 2 , and it is Schur stable (all roots in the open unit disk) if and only if 



1 1 



Pk < 0, k — 1, 2, . . . , 4 and p 6 = -P2P3P4 + P2P5 + < where 



Pi " 




' -1 


1 


-1 


1 


-1 " 




" 9o ' 


P2 




4 


-2 





2 


-4 




9i 


P3 




-6 





2 





-6 




92 


Pi 




4 


2 





-2 


-4 
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P5 




-1 


-1 


-1 


-1 


-1 
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The affine inequalities p&(:ci,X2) < 0, k = 1, 2, . . . , 5 define a polytope in the controller 
parameter plane (xx,X2) G M 2 , and the inequality ps(xi,X2) < defines a cubic region. 

In the case a = 0, with the following Gloptipoly 3 implementation of Steps 1-3 of Algo- 
rithm 1: 

mpol x y 2 

pi = -x(l)+x(2); 

p2 = -6*x(l)+4*x(2) ; 

p3 = -10*x(2)-4; 

p4 = -8+4*x(2)+6*x(l) ; 

p5 = -4+x(l)+x(2) ; 

p6 = 6*x(l)"2*x(2)+3*x(l)"2-10*x(l)*x(2)-2*x(l)-3*x(2)"3+6*x(2)~2+x(2) ; 
g6 = diff(p6,x); % gradient 
H6 = diff(g6,x); % Hessian 
% LMI relaxation order 

order = input (' LMI relaxation order = '); 
% build LMI relaxation 

P = msdp(min(y'*H6*y) , p6==0, pl<=0, p2<=0, p3<=0, p4<=0, p5<=0, ... 

g6*y==0, y'*y==l, order); 
°/„ solve LMI relaxation 
[status, obj] = msol(P) 

we obtain a negative lower bound obj = -3 . 5583 at the 2nd LMI relaxation, which is 
inconclusive. At the 3rd LMI relaxation, we obtain a positive lower bound obj = . 8973 
which certifies convexity of the stability region, see Figure [6j Since the stability region 
S = {x G M 2 : Pk{x) < 0, k = 1, 2, . . . , 6} is convex, we can optimize over it with standard 
techniques of convex optimization. More specifically, a recent result in [TT] indicates that 
any limit point of any sequence of admissible stationary points of the logarithmic barrier 
function f(x) = — Ylt=i \°&Pk( x ) is a Karush-Kuhn- Tucker point satisfying first order 
optimality condition. In particular, the gradient of f(x) vanishes at the analytic center 
of the set. Using Maple (or a numerical local optimization method) we can readily obtain 
the analytic center 2* w 0.57975, x* 2 ~ 0.13657 (five-digit approximations of algebraic 
coefficients of degree 17) corresponding to a controller well inside the stability region. 
Such a controller can be considered as non-fragile, in the sense that some uncertainty on 
its coefficients will not threaten closed-loop stability. 

Now for the choice a = —3/4 we carry on again our study of convexity of the stability 
region with the help of a similar GloptiPoly script. At the 2nd LMI relaxation we obtain 
a negative lower bound obj = -385 . 14 which is inconclusive. At the 3rd LMI relaxation, 
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Figure 6: Convex stability region (dark gray), with analytic center (cross) corresponding 
to a fixed-order controller. 

we obtain a negative lower bound obj = -380 . 88 which is also inconclusive. Eventually, 
at the 4th LMI relaxation, we obtain a negative lower bound obj = -380.87 which is 
certified to be the global minimum with status = 1. The point x 1 at which the minimum 
curvature is achieved is a vertex of the stability region, and the tangent at this point of 
the nonconvex part of the boundary is used to generate a valid inner approximation S, see 
Figure [3 Any point chosen in this triangular region corresponds to a stabilizing controller. 
We see that here the choice of the point of minimum curvature is not optimal in terms of 
maximizing the surface of S. A point chosen elsewhere along the negatively curved part 
of the boundary would be likely to generate a larger convex inner approximation. 

5.3 Optimal control with semialgebraic constraints 

In Model Predictive Control (MPC), an optimal control problem is solved recursively. 
This resolution is usually based on direct methods that consist of deriving a nonlinear 
program from the optimal control problem by discretization of the dynamics and the path 
constraints. Since the embedded software has strict specification on algorithm complexity 
and realtime computation, convexity of the program is a key feature [17J. Indeed, in this 
context, our convex inner approximation of the admissible space become valuable to speed 
up the computation even at the price of some conservatism. 

In open-loop control design, convexity of the problem is a matter of concern especially 
when the optimal control problem is part of an MPC procedure. In this case, the optimal 
control problem is solved mostly using direct methods that transfom it into a parametric 
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Figure 7: Convex inner approximation (dark gray) of nonconvex fourth-order discrete- 
time stability region (light gray). 

optimization problem. Convexity permits to limits the complexity of the resolution and so 
reduces the computation time of an optimal solution. Unfortunately, in dynamic inversion 
techniques based on differential flatness, the generally convex constraints on the states 
and inputs are replaced by nonconvex admissible sets in the flat output space, see [17] 
and reference therein for details. Thus, in such a method, it is necessary to design inner 
convex approximation of the admissible subset to develop a tractable algorithm [TS] . 

Consider the following optimal control problem 

min^ 
s.t. 



The objective of this problem is to steer the linear system from an initial state to a final 
state in a fixed time inside the admissible state subset S defined e.g. by the waterdrop 
quartic defined in section 14.31 

S = {i6R 2 : px(x) = x\ + x\ + x\ + x\< 0}. (6) 

We describe thereafter a classical methodology for solving the previous optimal control 
problem using flatness-based dynamic inversion, see [131 HH US] for other examples. As 
the dynamics are linear and fully actuated, dynamic inversion can be used to develop 



u 2 {t)dt 
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Pi(x 
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= x 
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Xf 



< 0. 
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an efficient algorithm for the considered problem [19] . Thus, the system trajectory can 
be parametrized by a user-specified sufficiently smooth function x\{t) = f(t) so that 
X2{t) = f(t) and u(t) = fit). The function f(t) is classically described by a chosen basis 
b(t) and the associated vector of weighting coefficient a such that 

f(t) = ^a k b k (t). 

k 

In order to derive a finite dimensional program, the admissible set constraint is discretized 
and enforced at a finite number of time instants {ti}i=i,...,jv such that t < ti < t 2 < ■ ■ ■ < 
tN < tf. Since S is nonconvex, we obtain a finite-dimensional nonlinear nonconvex 
program: 

min a J2i u2 ( a ,ti) 

s.t. x(a, to) — Xq, x(a,tf) — Xf 

Pi{x{a,ti)) <0, i = l,...,N. 

The inner approximation S calculated previously in section 14.31 is given by S = {x G R 2 : 
Pi(x) < 0, Pi{x) = gi(x 1 )(x — x 1 ) < 0, pz{x) = gi(x 2 )(x — x 2 ) < 0} where gii^x 1 ) and 
gi(x 2 ) is the gradient of p\(x) evaluated at x = x 1 and x = x 2 , respectively. The use of 
the inner approximation S as admissible subset leads to the following convex program: 

min a Y.i u2 ( a ^i) 

s.t. x(a,t ) = x , x(a,tf)—Xf 

p 1 (x(a,t i )) < 0, p 2 (x(a,ti)) < 0, pz{x{a, U)) < 0, i = l,...,N. 

In the following, we set t = 0, x = [0.3000, -0.8000] and t f = 2.5, x f = [-0.3000, -0.8000]. 
The time function f(t) is a 5-segment-piecewise polynomial of the 4th order (degree 3) 
defined on a B-spline basis. We run both programs for different values of N. In tabled] we 
compare the computation times and optimal costs. See Figure [8] for the state trajectories. 
For this example, we observe the positive effect that convexity has on the reduction of 





N 


10 20 50 100 200 500 1000 


CPU time [s] 


Convex 


0.029 0.045 0.056 0.058 0.102 0.225 0.498 




Nonconvex 


0.109 0.195 0.199 0.409 0.513 0.836 1.46 


Optimal cost 


Convex 


1.65 1.67 1.68 1.69 1.68 1.68 1.68 




Nonconvex 


1.52 1.52 1.52 1.52 1.52 1.52 1.52 



Table 1: Computation times and optimal costs of the nonconvex and convexified optimal 
control problems, as functions of the number N of discretization points. 

the computational burden, balanced by the relatively small loss of performance. 

6 Conclusion 

We have presented a general-purpose computational algorithm to generate a convex inner 
approximation of a given basic semialgebraic set. The inner approximation is not guar- 
anteed to be of maximum volume, but the algorithm has the favorable features of leaving 

15 



Figure 8: Optimal trajectories (bold) in nonconvex admissible set (left) and in convex 
inner approximation (right). 

invariant a convex set, and preserving convex boundaries while removing nonconvex re- 
gions by enforcing linear constraints at points of minimum curvature. Even though our 
initial motivation was to construct convex inner approximations of stability regions for 
fixed-order controller design, our algorithm can be used on its own for checking convexity 
of semialgebraic sets. 

Each step of the algorithm consists in solving a potentially nonconvex polynomial opti- 
mization problem with the help of a hierarchy of convex LMI relaxations. For this we 
use Gloptipoly 3, unfortunately with no guarantee of a priori computational burden, even 
though in practice it is observed that global optimality is ensured at a moderate cost, as 
soon as the dimension of the ambient space is small. Numerical experiments indicate that 
the approach may be practical for ambient dimensions up to 4 or 5. For larger problems, 
we can rely on more sophisticated nonlinear or global optimization codes [TS], even though 
this possibility has not been investigated in this paper. Indeed, our main driving force is 
to contribute with a readily available Matlab implementation. 

Our algorithm returns a sequence of polynomials such that the intersection of their sub- 
level sets is geometrically convex. However, the individual polynomials (of degree two 
or more) are not necessarily convex functions. One may therefore question the relevance 
of applying a relatively complex algorithm to obtain a convex inner approximation in 
the form of a list of defining polynomials which are not necessary individually convex. 
A recent result of [IT] indicates however that any local optimization method based on 
standard first-order optimality conditions for logarithmic barrier functions will generate 
a sequence of iterates converging to the global minimum of a convex function over convex 
sets. In other words, geometric convexity seems to be more important that convexity of 
the individual defining polynomials. 

Indeed, if convexity of the inner approximation is guaranteed in the presented work, con- 
vexity of the defining polynomials would allow the use of constant multipliers to certificate 
optimality in a nonlinear optimization framework. Instead, with no guarantee of convex- 
ity of the defining polynomials, the geometric proprety of convexity of the sets is more 
delicate to exploit efficiently by optimization algorithms. 
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Finally, let us emphasized that it is conjectured that all convex semialgebraic sets are 
semidefinite representable in [3], see also [TO]. It may then become possible to fully 
exploit the geometric convexity of our inner convex through an explicit representation 
as a projection of an affine section of the semidefinite cone. For example, in our target 
application domain, this would allow to use semidefinite programming to find a suboptimal 
stabilizing fixed-order controller. 
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